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Abstract 

The propagating chemical fronts found in cubic autocatalytic reaction- 
diffusion processes are studied. Simulations of the reaction-diffusion 
equation near to and far from the onset of the front instability are 
performed and the structure and dynamics of chemical fronts are stud- 
ied. Qualitatively different front dynamics are observed in these two 
regimes. Close to onset the front dynamics can be characterized by a 
single length scale and described by the Kuramoto-Sivashinsky equa- 
tion. Far from onset the front dynamics exhibits two characteristic 
lengths and cannot be modeled by this amplitude equation. An am- 
plitude equation is proposed for this biscale chaos. The reduction of 
the cubic autocatalysis reaction-diffusion equation to the Kuramoto- 
Sivashinsky equation is explicitly carried out. The critical diffusion 
ratio 6c, where the planar front loses its stability to transverse pertur- 
bations, is determined and found to be 5c = 2.300. 
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I. INTRODUCTION 



Propagating fronts separating regions with different characteristics occur in many 
physical contexts. Often such fronts have a comphcated structure of their own as 
in the fractal forms that arise in some diffusion-limited aggregration processes or 
viscous fingering [||]. In this article we consider chemical fronts separating regions 
of distinctly different chemical concentrations and study the nature of the spatio- 
temporal dynamics that these fronts exhibit. In these systems, under appropriate 
conditions, a planar front may not be stable and complex, even chaotic, front dynam- 
ics can arise. For example, this is the situation often encountered in propagating 
flame fronts which have been studied extensively f^. The complex dynamics of 
fronts can also underlie and determine the character of spatio-temporal chaos and 
dynamics in a variety of other chemical contexts 0|. 

We study a specific system here, the fronts in the cubic autocatalysis reaction 
A + 2B —* 3B, but several parts of our analysis are general and should apply in 
other situations. These cubic autocatalysis fronts have many features in common 
with propagating flames 0. An experimental realization of such chemical front in- 
stabilities occurs in the iodate-arsenous acid reaction carried out in a gelled medium 
to suppress fluid flow 0. The experiments show patterned fronts, much like those 
in quadratic and cubic mixed-order models. 

For the cubic autocatalysis model, when the diffusion coefficient of the fuel A 
is sufficiently larger than that of the autocatalyst B, the planar front is unstable 
to transverse perturbations. We have investigated the nature of the resulting front 
dynamics as a function of the diffusion ratio Da/Db in large systems. Several 
results have emerged from our study of this system. For small enough diffusion 
ratios within the unstable regime the front exhibits complicated dynamics with 
statistically stationary properties and is characterized by a single length scale. Far 
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beyond the instability point a new chaotic regime is encountered, termed biscale 
chaos, in which there are two characteristic length scales. 

We have carried out a detailed analysis of the front dynamics in terms of am- 
plitude equations. Just beyond the instability threshold the front dynamics is de- 
scribed in terms of the Kuramoto-Sivashinsky equation We document this 
relation through a statistical characterization of the interface. The regime far from 
the instability point where biscale chaos is observed cannot be described in terms 
of the Kuramoto-Sivashinsky equation and we construct an amplitude equation in 
which the nonlinear mode coupling is generalized. This amplitude equation is able 
to capture the principal qualitative features of the biscale chaos. This part of the 
analysis may generalize to other situations. 

We also give a detailed reduction of the cubic autocatalysis reaction-diffusion 
equation to the Kuramoto-Sivashinsky equation. The coefficients that enter in this 
equation are estimated analytically using solutions for the right and left eigenvectors 
of the eigenvalue problem that enters the analysis, and a numerical scheme is devised 
that allows the computation of these coefficients for any values of the diffusion 
coefficients. As a by-product of this analysis we may easily determine the critical 
diffusion ratio where the planar front becomes unstable. 



II. CUBIC AUTOCATALYSIS FRONTS 

Consider the autocatalytic reaction A + 2B — > 3B described by the reaction- 
diffusion equation, 

= -af3^ + DAAa , 

ot 

^ = a/32 + ^^^^ ^ 
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where a{t,r) and (3{t,r) are the (scaled) concentrations of the A and B species, 
respectively, with diffusion coefficients and Db- In one dimension, with suitably 
defined initial conditions ||^, this system supports isothermal chemical fronts 0. 

In two space dimensions the planar front is unstable to transverse perturbations 
when Da is sufficiently larger than Db- The origin of such instabilities and the 
nature of the resulting front dynamics have been the subject of earlier studies [Q. 
Horvath et al. |T^ performed simulations of the cubic autocatalysis model (|1]) and 
investigated the bifurcation structure as a function of the system length. For small 
system sizes, when the front possesses a single minimum, the onset of chaotic be- 
haviour was observed to occur through a period doubling cascade in the temporal 
dynamics of the minimum. For large system sizes the resulting front dynamics is 
best analysed in statistical terms. We now briefly describe and quantitatively char- 
acterize the nature of the front dynamics for large system lengths as a function of 
the diffusion coefficient ratio, 6 = D^/ Db- 

We consider a two-dimensional system which is infinite in the x direction and 
has length L along where periodic boundary conditions are applied. The initial 
conditions are taken to be {a(0,x, ?/) = 1, /?(0,x, = | x < —w/2,Wy} and 
{a{0,x,y) = 0, P{0,x,y) = I \ x > w/2,\/y}- The region {x,y \ — w/2 < x < 
w/2,0 < y < L} was divided into segments of length i along y and each segment was 
assigned the values (a(0, x, y), P{0, x, y)) = (1, 0) or (0, 1) with probability p = 1/2. 
It is convenient to represent the front dynamics in terms of the isoconcentration 
profiles, hA(t,y) and hB(t,y) defined by a(t,hAit,y),y) = ar and P{t,hB(t,y),y) = 
Pr, where and (3r are reference concentrations. Henceforth we focus on hB{t,y) 
measured relative to its average value at time t and denote this quantity by h{t, y). 
We note that for these initial conditions this average value moves with the minimum 
front speed, which is selected from the continuum of allowed front speed values. 
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Figure |I|(a) shows a space-time plot of the minima of h{t,y) for 6 = 5. In the 
integration |Tl| of (0) we have scaled the length variables so that Da = 1 and 
Db = 6^^. After an initial transient period during which the interface develops, a 
moving front with statistically stationary properties is formed. In this stationary 
regime the spatio-temporal behaviour can be described in terms of the dynamics of 
the minima which play the role of "particles" in the system. The particles collide and 
coalesce and new particles are born so that the average density of particles per unit 
length remains constant. The front dynamics may be characterized quantitatively 
by the space-time correlation function 

Cit, y) = (L-i r dy'[h{t, y') - h{t, y + y')f) , (2) 

JO 

where the angle brackets signify an average over realizations of the front evolution 
starting from the random initial conditions given above. We also consider the time 
average of C{t, y) in the stationary regime, 

C{y) = T-' f"''^dt'C{t',y), (3) 

where to is a time longer than the transient period. The correlation function 
is sketched in Fig. ^ for several values of the time t and shows the development of 
correlations during the transient period and the approach to the stationary regime. 
The function C{t,y) exhibits short range correlations with a tendency to reach a 
nearly parabolic form, seen in the C{y) curve, which is characteristic of diffusive 
transverse front dynamics and will be discussed further below. The persistent short 
range correlations are due to the existence of a characteristic length i* arising from 
the front instability. 

This type of behaviour persists up to about 5 = 6 when a new type of chaotic 
front dynamics, distinguished by the appearance of a second length scale, is observed. 
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The second length scale is easily discerned in both the space-time plot of the minima 
of h{t, y) shown in Fig. |I](b) for (5 = 8 and the power spectrum, 



The power spectrum is presented in Fig. ^ for both 5 = 5 and 5 = 8. For 5 = 5 
one sees a single peak corresponding to the wavenumber of the most unstable mode. 
However, for 5 = 8 the minimum in the power spectrum has filled in, indicating the 
growth of modes with smaller wavenumbers and thus the appearance of structure 
on longer length scales. We now turn to an analysis of these results in terms of 
amplitude equations and provide a foundation for the phenomenological picture of 
the front dynamics. 



A description of the origin of the front dynamics is most conveniently given in 
terms of amplitude equations for the front profile derived from the reaction-diffusion 
equation (|lD. We analyse the dynamics close to the instability point in terms of the 
Kuramoto-Sivashinsky equation and show how this equation must be generalized in 
order to explain the regime far beyond the instability where a second characteristic 
length appears. 



We present an adaptation of Kuramoto's derivation for the case of small 
amplitude perturbations which yields the amplitude equation that will be used in 
the subsequent analyses. We shall confine our attention to the question of how the 
presence of a second dimension affects the dynamics of the general reaction-diffusion 

system. 




to 



(4) 



III. ANALYSIS OF FRONT DYNAMICS 



A. Dynamics of small perturbations 
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^ = F(z) + DAz . 

at 



(5) 



Here z is a vector of concentration fields and F(z) is a vector-valued function de- 
scribing chemical reactions. We assume that in one dimension (1) possesses a stable 
solution with a propagating front profile z(t, x, y) = zo(x — ct), where c is the veloc- 
ity of the front. Following the treatment by Kuramoto |jl2| we seek the dynamics of 
a perturbed solution in the form, 



z(^, X, y) = zo(^ + Mt, y)) + <Pi{t, 



(6) 



i>0 



where ^ = x — ct. As is customary, we use an abstract notation (^|ui) = Ui(^). The 
vectors |uj) are the solutions of the eigenvalue problem. 



d'^ d 
OT(zo) + + c— . 



\ui) = C\ui) = Xi\ui) 



(7) 



where PF is the Jacobian of the vector-function F and C is the operator in square 
brackets. Here a hat signifies an abstract operator while the corresponding quantity 
without a hat is its ^ representation. We expand F(z) in a Taylor series near zq 
and neglect terms of second order and higher. Since the main features of the front 
dynamics are embodied in T>F, the truncation of the Taylor series does not discard 
any important information. Due to the translational symmetry of equation (j^) we 
have a straightforward solution of (0), namely |uo) = (9|zo)/9^ corresponding to a 
zero eigenvalue. Taking this solution into account we write 



d 



^|ui) = c|uo)+F(zo) + D-^zo + (V0o) D-^|uo) + 



dt 



i>0 



W|u,) + D-^|u,)+c-^|u,) 



(8) 



Using the identity — c|uo) = F(zo) + D — ^|zo) and definition (0) we arrive at the 



expression 
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.2^ d 



-^\ui) = A,0,|u,) + A0,D|u,) + (V0o) D-.|uo) . (9) 

Here and below we use the Einstein summation convention. Alternatively, after 
multiplication by (uj|, (H) may be represented by a set of coupled equations 



at ...... 



Xi^i + (ui|D|u,)A0, + (ui|D — |uo) (V0o)' • (10) 



This set of equations can be formally solved for modes (pi {i ^ 0) by treating 
00 as an independent function and applying Duhamel's principle to the resulting 
system of inhomogeneous linear equations. We find 



t 

4>{t) = e^*0(O) + j dse^^'~'^ [aA0o(s) + b (V0o(s))' 
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In (|TT]) we use the following notation: the matrix operator W has elements 

Wij = 5ijXi + (ui|D|uj)A, the vectors a and b have elements = (ui|D|uo), 

d 

hi = (uj|D — ^|uo), respectively, and is a vector with elements 0j (i > 0). 

Assuming that exp(Wt) decays rapidly compared to 0o we may perform the 
integration in ( pT]) and substitute the result into the equation for (pQ to obtain 

ao - E c,{W-%jajA\ A0o + [bo - E c.{W-i},,6,A | (V0o)' , (12) 

i,j>0 J \ i,j>0 J 



dt 



where c, = (uo|D|uj). This is the generalized amplitude equation which will form 
the basis of the analysis of the biscale front chaos given below. 

If the dynamics is described by small k modes one can take W^^ to be a diagonal 
matrix 6ij/Xi and by omitting second order terms in the coefficient of (V0)^ one 
obtains IT! 



^0 / IT^I \ A J, , / IT^ ^ I \ /V7^ n2 V- (uo|D|ui)(ui|D|uo) . 2 



uo|D|uo)A0o + (uo|D-.|uo) (V0o)' - E ' ' ' V ^' ' ' A ■ (13) 



dt ^ - ' - ^ - A 



The coefficient of the (V0o)^ term is obtained from the following identity: 
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= 



(uo|i:|0(e|uo)-(uo|0(^l^|uo) 



2(uo|0(e|D^|uo)+ (^c-D^^ (uolO(eiuo) 



(14) 



The above formula vanishes for all x so that after integrating over the entire do- 

c 

main we obtain (uo|D^|uo) = — - and observe that the expression has a universal 
character. This gives the Kuramoto-Sivashinsky equation, 



dt 



I (V0o)' 



(15) 



with 



z/ = (uo|D|uo) 



(16) 



and 



E 

j>0 



(uo|D|ui)(uj|D|uo) 

A,; 



(17) 



A stability analysis of (0) shows that the planar front is unstable with respect to 
long-scale, small-amplitude perturbations when z/ < 0. 

In the next subsection we analyse the cubic autocatalysis front dynamics in terms 
of both the generalized amplitude equation (|12D as well as the Kuramoto-Sivashinsky 
equation (|T5|). 



B. Amplitude equation description of front dynamics 

1. Dynamics near instability 

The Kuramoto-Sivashinsky equation provides insight into the nature of the cubic 
autocatalysis fronts close to the instability point, and we briefly comment on this 
connection. Simulations of the Kuramoto-Sivashinsky equation produce height cor- 
relation functions and space-time plots (cf. Fig. ^) of the front dynamics similar to 
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those of the cubic autocatalysis model. In terms of the Kuramoto-Sivahinsky equa- 
tion, one sees that the front dynamics for 5 > 5c is determined by the instabihty at 
long wavelengths, due to the fact that u < in the A0o term, and the dissipation 
at short wavelengths controlled by the A^0o term with positive k. The dispersion 
relation of the Kuramoto-Sivashinsky equation is u!{k) = —uk^ — nk"^ with maximum 
at km = (— z//2/t)^/^. The long wavelength structure drives the dynamics of the ex- 
trema ( "particles" ) whose collision dynamics provides the dissipation of the energy. 
The characteristic distance between extrema can be related to the wavevector by 
i* = 2TT/km- Using a zeroth order approximation to the u and k coefficients given in 
the next section, for 5 = 5 we obtain i* = 30 which is is comparable to the value of 



i* = 39 found in the numerical simulations. Using the methods given in Sec. |^ it 
is possible to obtain u and k accurately and improve this estimate. For the present 
purposes this rough estimate plus the appearance of the space-time plot ||14| in Fig. ^ 



suffice to confirm the general character of the Kuramoto-Sivashinsky mechanism for 
the chaotic cubic autocatalysis front dyanmics for large system lengths. 

It has been argued |]T5| that the long-scale dynamics of the one-dimesional 



Kuramoto-Sivashinsky equation can be described by the Kardar-Parisi-Zhang equa- 
tion p!6[] . 



ahMl = v^L^^i^vMt,v)? + at,v), (18) 

where c is the front speed, D is a diffusion coefficient and ^(t) is a Gaussian white 
noise source with correlation function {C,{t,y)C^{t',y')) = 2T6(t — t')6{y — y'). Al- 



though there is debate about some aspects of the reduction in higher dimensions [17 



numerical simulations |jT8| of the scaling properties of the Kuramoto-Sivashinsky 
equation have confirmed this reduction in one dimension. The diffusive character 
of the cubic autocatalysis front dynamics discussed in Sec. |I| is consistent with 
such a description and provides additional confirmation of the applicability of the 
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Kuramoto-Sivashinsky equation close to the onset of instability. 



2. Biscale chaos far from instability 

Far from instability, for 6 ^ 6c, we have seen that the cubic autocatalysis model 
presents new features. The space-time plot (Figure 0(b)) shows that besides the rela- 
tively short-lived trajectories, which are similar to those of the Kuramoto-Sivahinsky 
equation, there exist long-lived, stable patterns. We have verified that this type of 
biscale chaotic front dynamics, which involves a second characteristic length scale, 
cannot be described by the Kuramoto-Sivashinsky equation. One may seek a gener- 
alization of the Kuramoto-Sivashinsky equation in which the coefficients of the A(f)o 
and A^0o terms are modified; this in effect changes the dispersion relation and if it 
has the form sketched in Fig. ^(a) with two extrema then one might observe biscale 
spatio-temporal dynamics. We note, however, that direct numerical calculation of 
the fastest growth mode and the growth rate corresponding to it (cf. Appendix for 
the method used to carry out this calculation) shows that the dispersion relation 
remains similar to that of the Kuramoto-Sivashinsky equation even for 6 ^ 6c (see 
Fig. ^ for a schematic representation of the situation). Thus one must seek a de- 
scription of this new type of front dynamics outside the usual Kuramoto-Sivashinsky 
equation or modifications of it that involve a generalized dispersion relation. 

We now describe a model that captures the qualitative features of this phe- 
nomenon and is consistent with the amplitude equation (|I^) derived above. Since 
the origin of the new dynamics does not reside in terms that affect the dispersion 
relation, one is led to examine the nonlinear term in (0) which, for future reference, 
we rewrite in terms of its Fourier transform: 

dMt, k) ^^^^^^(^^ A;) + ^ 5: k'{k - k')Mt, k')Mt, k-k'), (19) 



dt ' ^"'^^ ' ^ ' 47r 



11 



where g{k) is given by 

- ^gik) = &o - E cdWik)-%b^e . (20) 

i,j>0 

Recall that 60 = — c/2. The dynamical system governed by d(j)/dt = (V*/))^ conserves 
the quantity J (VcpYd^ as can be seen easily from the identity ||T9| : 

j- 1 {Vct^fdi = -2 1 VV(V0)^rfe = -ll V(V0)^de , (21) 

so that below we shall refer to the linear and nonlinear terms in (p!9[) as dissipation 
and mixing terms, respectively. Multiplying (|l^) by (f)o{t, ~k) and averaging the 
resulting equation over time we have for the stationary power spectrum 

= uj{k)E{k) + ^ E ^'(^ - (<^o(t, A;')0o(t, k - k')Mt, -k)) . (22) 

From this equation it is clear that the stationary properties depend on the ratio 
uj{k)/g{k) so that for g{k) as in Fig. |^(b) we see that the growth of the mode with 
wavenumber k^ is strongly suppressed. While the above arguments do not constitute 
a rigorous "proof" they serve as a guide to the proper form of the mixing coefficient. 
Technically, the above scheme is accomplished in a natural way when a zero of W(fc) 
comes close to real axis in the vicinity of the maximum of the dispersion relation. 
In this case behaviour of g{k) is dominated by the pole of W^^(A;) and the mixing 
coefficient may be approximated by 

[k - koY + 7^ 

where the damping term can be rationalized in terms of the continuous spectrum of 

K 

In the simulation of this equation we scale space, time and 0o so that z/, k and 
c/2 are unity and take the following values for the parameters that appear in g{k): 
e = 0.2, fco = 0.6 and 7^ = 0.005. The results are presented in Fig. ^(b). In 
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the figure one sees the existence of the structures with two different length scales, 
the main signature of biscale chaos. The generalized amplitude equation does not 
capture all of the details of cubic autocatalysis front dynamics far beyond the onset 
of instability. For example, the simulations of cubic autocatalysis show a broader 
distribution of longer length scales and, as is quite evident in Fig. |1|, regions where 
there are no extrema. The full descriptions of these features may be specific to cubic 
autocatalysis and require a more detailed analysis. However, the main feature, the 
emergence of a second length scale, is described by the model. 

We note the similarities in transitions to ordinary and biscale chaotic regimes. In 
the former case the transition takes place when (uo|D|uo) passes through zero and 
in the latter we require a zero of {SijXi — {ui\Y)\uj)k'^} to cross the real axis. The 
correspondence between the pole of and the existence of structures with dif- 
ferent wavenumbers suggests that this phenomenon can arise in other contexts. The 
existence of more complicated structures, say with three or four different wavenum- 
bers, is unlikely in view of the continuity of the spectrum and the fact that the 
dispersion relation is a smooth function of k. Consequently, the amphtude equation 
(p!9|) could serve as the basis for the analysis of such phenomena in systems other 
than the cubic autocatalysis reaction. 

IV. PARAMETERS IN KURAMOTO-SIVASHINSKY EQUATION 

Thus far we have not related the parameters in the Kuramoto-Sivashinsky equa- 
tion to those in the cubic autocatalysis model (|l]) and in this section we determine 
u and K. We first present an approximate analytical calculation of these parameters 
based on the exact solutions of the right and left eigenvalue problems for |uo)and 
(uo|, respectively, for case of equal diffusion coefficients pO|. This provides insight 
into the structure of the eigenvalue spectrum and the corresponding eigenvalues. 
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We then describe a general numerical scheme which can be applied to other models 
as well. In this section we choose a different space scaling and set = 1 + and 
Db = ^ — since this simplifies the computations. 



A. Case of equal diffusion coefficients 

When the diffusion coefficients are equal (/i = 0) the reaction-diffusion equation 
(HI) supports an additional relation involving the a and (3 concentrations due to 
conservation of the total mass at each point of space, a + (3 = 1, given that at the 
initial time this relation is satisfied at every space point. We can rewrite system (|l]) 
{Da = Db = 1) in the form 

^ = -ail-af + Aa. (24) 
at 

Local initial input of the autocatalyst induces the creation of a front, propagating 
with the minimal velocity 0, whose profile is 

ao(0 = (l + exp(-ce))-\ (25) 

where c = l/\/2 is the velocity of the front. The form of the [5 profile follows from 
the mass conservation relation. Henceforth, when confusion is unlikely to arise, we 
drop the subscript on a and j3 to simplify the notation. 

Next we discuss the spectrum of the eigenvalue problem (^ for the cubic auto- 
catalysis problem when the diffusion coefficients are equal. Using the explicit form 
of PF(zo) derived from (|1|) we have 

+ - /5']Xi - 2al3X2 = XX, , 

/^'^i + + + 2«/3]^2 = AX2 , (26) 

where we have written the components of a general right eigenvector u as u = 
(Xi,X2). We shall show that the spectrum is obtained from two distinct equations. 
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First, by adding the two lines of (p6D , we obtain a linear differential equation with 
constant coefficients: 



&^ d_ 



X = XX 



(27) 



where X = Xi + X2. After a change of variables X = exp(— c^/2)y we obtain 
9^/9^^ — Y = \Y. The eigenvalue problem for Y has a continuous spectrum 
00, — . In these considerations we have omitted the case Xi + X2 = 0, which 
satisfies the eigenproblem for X automatically. Using this relation we obtain the 
following equation for the other part of the spectrum: 

Xi = AXi . (28) 

After the same change of variables Xi = exp(— c^/2)Yi we get a "Schrodinger" 
equation with potential V = —2a[3 + (cf. Fig. ^(a)) and energy levels 
Ei = /A — Aj. We see that this equation has a continuous spectrum starting from 
zero as well as discrete spectrum. The WKB method may be used to obtain approx- 
imations for the eigenvalues and eigenfunctions. Applying the Bohr-Sommerfeld 
quantization rule we get Eq = —0.1066 for the first energy level, which is in good 
agreement with the theoretical value — = —0.125. This indicates that the WKB 



approximation is reliable and it predicts that the potential well is too shallow ||2T 



to have other discrete energy levels. So we conclude that the spectrum consists of a 
zero eigenvalue and doubly degenerate continuous spectrum lying in (^—00, — . 

The right eigenvector |uo) = (9|zo)/9^ was discussed in Sec. |I|, but the left eigen- 
vector (uqI is more difficult to obtain. We now construct the explicit form for (uo| 
for the case of equal diffusion. Letting (uo|^) = Uq(^) = {X^,X2) the eigenvalue 



problem for Uq(^) takes the form 



XI - (3'^ {XI - X*) 
X* - 2«/5(X* - X*) 








(29) 



15 



so that the difference X* = — X2 obeys 



This can be cast into the form of ( pSD after the change of variable X* = Aexp{cC,)Z, 
where A is some constant defined by the normalization condition (uq|uq) = 1: 



00 1 
—X*d^ = cA J c?dot 



cA 



(30) 



hence, we have 



X* = 3(l + exp(-cO)^ 



(31) 



Performing the change of variable \1/ = exp(— c.^/2)X^ in (|29D and using (^) we 
obtain 



92 



^ = 3/32exp(-ce/2)(l + exp(-ce))- 



(32) 



de ' V2 

The Green's function for the operator (9^/9^^ — is given by Q{i) = 
— exp(— c|^|/2). Convolution of Q with the right-hand-side of (|3^) yields a solu- 
tion for XI, 

(1 + exp(-cO)^ 
and the difference between X^ and X* gives Xg, 

4 + 2exp(-<^) 
^^'«*--(l+exp(-c«))^- 

These functions are plotted in Fig. P(b). The functions have limiting behaviour 
X*(+oo) = -1 and X*(+oo) = -4. 

System (p9D, for any diffusion coefficient ratio, has another zero left eigenfunction 
V* = (1,1), which appears due to conservation of the total mass of the reagents. 
The right eigenfunction V corresponding to V* is 
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3 + 3 exp(— G^) 
2 



2ce + 3 \ 
l + exp(cOy ' 
2c^ + 3 



(35) 
(36) 



3 + 3exp(— c^) Y l+exp(c^)y 

While both the "bra" vector corresponding to spatial invariance and the "ket" vector 
obtained above do not decay to zero at ,^ = oo their scalar product does and the 
product is given, after change of variable, by: 



1 

(V|uo) = V2j 



^o? In ■ 



a 



1 — a 



-2(l + 3a)(l - a) 



(ia = 



(37) 



Since the above integral vanishes these considerations insure us that we have an 
orthogonal system of "bra" and "ket" vectors. Moreover, the fact that one product 
is finite and the other is not removes an ambiguity in the definition of the left 
eigenfunction. We believe the same situation holds for arbitrary /x. 



1. Estimate of v for fJ- ^ 

For equal diffusion coefficients D is the identity matrix and it follows that u = 
(uo|D|uo) = 1. However, using |uo) and (uo| computed for /i = 0, we may obtain 
a first order approximation for u for unequal diffusion coefficients. From (|16D we 
write: 



u = (uo|D|uo) = 1 + ;u(uo|I|uo) 



(38) 



where I 
get 



^ ^ 



. Using the explicit, equal-diffusion forms for (uo| and |uo) we 



— {XI + X*) d^ = - J («' + 4«) da = -- 

-co 

and z/ = 1 — 7/i/3. 



(39) 
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2. Estimate of k for /i 7^ 



When the diffusion coefficients are equal and D is the identity matrix we have 
K = in view of the orthogonahty of the eigenfunctions. However, as in the preceding 
subsection we can obtain an estimate for k, for arbitrary /i using the solution of the 
eigenvalue problem for fi = 0. To obtain this estimate we assert completeness of the 
set of eigenfunctions; i.e., ^ |uj)(uj| = 1. Using this assumption we obtain 



^ (uo|A|ui)(u,|5|uo) / \ / iRi \/ I \ 

2^ ^ = (uo|fi|x) - (uo|fi|uo)(uo|x) 



A, 



(40) 



where |x) is a solution of the equation 

|x) = A\uo) - (uo|A|uo)|uo) . 



(41) 



In this gereral expression we let A = B = JD and use (uo|D|uj) (uj|D|uo) = 
yU^(uo|I|uj)(uj|I|uo) in (p!]). Using A = D and adding the two lines of pl|) we 
have 



d_ 



{xi + X2) = 2a^ 



(42) 



with = da/dC,- Solving the above equation with appropriate boundary conditions 
we get 

21og(l + e'=«) 



Xi + X2 = — 



(43) 



c 

Equation (|4ll) can be solved for xi to obtain (we use in the next equation the exphcit 
form for (uo|I|uo)) 



xi = 2aP{xi + X2) + 



(44) 



This equation takes an especially simple form after a change of variable ^ — a. 



The derivative possesses the form 



d d 

— = ca{l — a-iid (|4^) after simplifications 



is tantamount to 



, , ^ d / 8/5 (1 -a)ln(l-a 



da^ da ^ ^ c \\2 a J 

(45) 

We know that af3 is a solution of the homogeneous equation. By the substitution 



Xi = aPx 'we decrease the order of (|4^). Solving the resulting linear first order 
equation we obtain the solution 

xi = — ( (i _ a^) Mn(l - a) + a + — ) + 7a (1 - a) ln(^^) ) . (46) 

From the form of the solution it is clear that the solution decays faster then e"'^'^'^', 
< 7 < 1, and thus the scalar product is well defined. To integrate we use the 
identity 

Qd^ = -j , (47) 

c J ail — a) 

and integration gives 

(uo|D|ui)(ui|D|uo) 264 + 407r2 ^ 
^ = 2^ T = — 7^ — ^^ ■ (48) 

Finally, collecting the above information, we arrive at the following expression 

for the Kuramoto-Sivahinsky equation with approximate coefficients deduced from 

the cubic autocatalysis model @): 

d(f)o ,^ 7 , 264 + 407r2 c ,^ , ,3 /.^n 

^ = (1 - -/i)VVo ^ /u'VVo - 2(V</'o)' • (49) 



B. Numerical computation of 1/ for 7^ 

In order to determine u exactly, within the present approach, one must construct 
a scheme to determine the left eigenvector Uq, the solution of C^u* = Au*, for A = 
and arbitrary /i. For the cubic autocatalysis problem the left eigenvectors are given 
by the solutions of 
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(9^ (9 



^1 + /3^^2 = A^i 

^2 = 



(50) 



where we have written a general left eigenvector as u* = (\l/i, \l/2)- The solution of 

this eigenvalue problem in terms of a complete basis entails the solution of a difficult 

sparse, non-symmetric, non-tridiagonal matrix eigenvalue problem. In order to avoid 

solving the eigenvalue problem by matrix methods we use the following device. We 

consider an auxiliary problem that originates from the system of equations (^0|). In 

d 



(M) we replace terms c— by the time derivative and take a and (3 to be in a moving 
frame with velocity — c (i.e., a{t, x) = ao{x + ct) , /3(t, x) = Po{x + ct) ): 



dt 
dt 



:i+/i)A^i-/32(^,-^2) 
:i - /i)A^i - 2«/?(^i - ^2) • 



(51) 



Eigenvalues of the adjoint eigenvalue problem C^u* = Au* coincide with those of (|^ 
and for nonzero eigenvalues their real parts are negative and separated from zero. 
Using this observation we deduce that all but one of the components exponentially 
decay so that the persistent solution of the system of partial differential equations 
(pip is proportional to (uo|. It is worth noting here that in order to obtain (uqI one 



should take unequal initial distributions for \I'i and \I'2 in (^TJ). 

Numerical computations of both (|l|) and ( pi]) were carried out using the following 

hybrid of the Runge-Kutta and Crank-Nicholson methods: 

z(t) z' + z(t) , 
^ ^ DA ^ + F(z(t)) , 

(52) 



z(t + r) - z(t) 



DA 



z(t + r) +z{t) 



z' + z(t) 



T 2 V 2 ; 

We use a non-uniform distribution of mesh points with density x/ (1 — x^)^^, where 
X is in the range [—1, 1]. The matrix coefficients and vectors were computed using 
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a finite element method with functions approximated by a sum of hat functions. 
During the simulations the profiles were monitored and constantly shifted to the 
origin. To do this the profiles were interpolated using a monotonicity-preserving, 
piecewise-cubic, Hermite interpolant to a set of data points For the velocity 
we use the following relation: 



oo 

J ap^d^. (53) 



Due to the fast decay of the rate of the reaction outside the reaction zone the infinite 
range of the integration can be replaced by the limit mesh points and the value of 
the definite integral of the interpolant is easily computed. 

The results obtained in simulations for the case of equal diffusion coefficients 
were compared with the analytical results from the previous section. The absolute 
errors for the velocity and for (uo|I|uo) are about 10~^ and 10"*^, respectively. This 
confirms that the numerical procedure is accurate. 

V. CRITICAL DIFFUSION RATIO 

While it is clear that the planar front is unstable for sufficiently large values 
of the diffusion ratio 6, no estimate was provided the critical diffusion ratio, 6c, 
where this instability takes place. There are a number of estimates of 5c in the 



literature. Horvath et al. [|T0| estimated 5c ~ 2.9 from direct simulation of the 
reaction-diffusion equation and found 5c = 2 from an application of Sivashinsky's 
method, which assumes an infinitely thin reaction zone, to the cubic autocatalysis 
problem. Zhang and Falle computed the dispersion relation directly and found 
5c ~ 2.0. Finally, Milton and Scott obtained 5c ~ 2.366 using the eikonal theory 
and 5c ~ 2.25 from numerical simulation. This quantity is difficult to estimate by 
simulation of the reaction-diffusion equation and its direct calculation presents a 
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number of interesting features. 



In the context of the amphtude equations discussed in Sec. 5c can be de- 
termined from the condition v = (uo|D|uo) = 0. This, in turn, requires the com- 



putation of the left and right eigenvectors of the eigenvalue problem The 
analytical estimate of u discussed above immediately gives an estimate of the criti- 
cal diffusion coefficient ratio. The critical value of /i, /Xc, is given by the solution of 
u = = 1 — Tjic/'i or He = 3/7. Since 6 = {1 + fi) / {1 — fi) we have 6c = 2.5. 

To determine the critical diffusion ratio exactly we use the numerical scheme 
outlined above to determine the left and right eigenvectors for arbitrary /i. In 
this way the critical diffusion ratio can be obtained iteratively. At each step the 
ratio from the previous run was taken and new profiles and the ratio — (uo|I|uo)^^ 
corresponding to them were determined. Iterations rapidly converge to the value 
He = 0.39397, which corresponds to diffusion coefficient ratio 6c = 2.300. 

We have also computed the dispersion relation to provide additional confirmation 
that 6c > 2. The details of the computation are given in the Appendix. The 
method we use allows us to put several thousand mesh points in the reaction zone 
(compared to about 40 in |2^). In the small wavenumber range the dispersion 
relation should have a dependence on k given by that of Kuramoto-Sivashinsky 
equation. The maximum growth rate for the dispersion relation uj{k) = —uk'^ — nk'^ 
is i/^/4k. We see that when 6 approaches 6c, so that z/ ~ ((5c — 6), the maximal 
growth rate rapidly decreases and can easily fall out of range of the precision of the 
calculation. Moreover, direct simulation shows the presence of a very slowly decaying 
tail of the autocatalyst and its length increases (cf. Appendix) when k approaches 
zero. So the calculations could also be affected by the limited system size. We 
also observe that least-square fittings of the data in Fig. 2 (a) of p4| gives u = .34 
and K = 1.7 indicating stability of the front. A further test of our computation 
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of the dispersion relation comes from the first order approximation of v in the 
Kuramoto-Sivahinsky equation when fi ^ fi^'- = fJ'c + ^A*- In this case we have 
u = 1 + yu(uo|I|uo) = —6fi/fic- The resuhs of our computations of the dispersion 
relation for the case /i = 0.45 {6 = 2.636) are given in Fig. 0. We see that near the 
origin the dispersion relation is well approximated by the curve u = 0.142/c^, where 
0.142 = Sfi/f^c. 

VI. SUMMARY AND DISCUSSION 

The chaotic dynamics of cubic autocatalysis fronts shows two distinct kinds of 
behaviour depending on the diffusion ratio Da/ Db- Close to the onset of the front 
instability the transverse structure can be described in terms of a single length scale. 
The characteristic length is associated with the wavelength of the most unstable 
mode and the dynamics can be viewed in terms of a particle picture. The front ex- 
trema can be associated with "particles" that collide and coalesce and new particles 
are created to maintain the average particle number. This dynamical regime was 
characterized statistically and shown to be described by the Kuramoto-Sivashinsky 
equation. 

For large values of Da/Db a new dynamical regime was found where the front 
is characterized by two length scales. We term this type of chaotic front dynam- 
ics biscale chaos. Biscale chaos cannot be described by the Kuramoto-Sivashinsky 
equation and we introduced an amplitude equation which is a generalization of this 
equation. This amplitude equation accounts for a modified nonlinear energy mixing 
among the modes. We identified this kind of behaviour with a pole near the real 
axis of a linear operator obtained in the course of the analysis of the dynamics of 
small perturbations. The observed phenomenon may prove to be quite general in 
light of the simple prerequisite for it, namely, the crossing of the real axis by a pole. 
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However, we believe that the coexistence of structures with three or more length 
scales is unlikely since it requires the operator to have more complicated properties. 
Questions about qualitative effects of this irregular dynamics on microscopic quan- 
tities (e.g., mean front velocity, diffusion coefficient of the interface, etc.) remain. 
Another point to be investigated is the possibility of describing the observed pat- 
terns in terms of the dynamics of interacting pulses, which correspond to small k 
modes, driven by noise. 

An explicit reduction of the cubic autocatalysis reaction-diffusion equation to 
the Kuramoto-Sivashinsky equation was carried out and the coefficients that enter 
in this equation were determined. Both approximate analytical values as well as a 
numerical scheme for their exact computation for arbitrary D^/Db were given. The 
numerical method copes with two problems that arise in such calculations, namely, 
the nonlocal form of the zero eigenvector and the degeneracy of the eigenproblem 
with eigenvalue equal to zero, which comes from the conservation laws for the system. 
The algorithm is easy to implement and proved to be reliable. Using these results 
the critical diffusion ratio where the front instability occurs was computed and found 
to be Da/Db = 2.300. 
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APPENDIX: NUMERICAL COMPUTATION OF THE DISPERSION 

RELATION 



The dispersion relation for the cubic autocatalysis model was computed in the 
following way: The perturbation of the front solution of (|^) may be written in the 
form 

z{t, X, y) = zo(0 + Xit, x) exp{tky) . (Al) 

We assume, when k is small, that the largest eigenvalue is separated from the others. 
Thus, after integration of the equation 

f = DFx + d(^-*^)x. (A2) 

for a sufficiently long time, only the most unstable mode survives. The system 



is the same as system (4.7) of p4| after proper scaling and change of variables 



daQ dhi, 



Xi = cii — -7—5 X2 = bi ;— in the notation of this paper [E^. After transients the 

ax ax 



solution of ([A2|) should propagate with the same velocity as the front and grow with 
a rate u!{k) and thus we can replace with u!{k)x — cx^- In @] the eigenvalue 



problem with % as a vector and uj{k) as eigenvalue was solved. The two approaches 
lead to the same result but the latter necessitates solving an eigenvalue problem 
for a non-symmetric, sparse matrix that is in general a difficult task. The direct 
integration if the system ( [A2D is more accurate and easier to implement. The results 
of the integration for fi = 0.45 are shown in Fig. 0. 

To decrease the integration time we used as initial conditions for x the eigen- 
vector obtained from the run for k — 6k and by employing numerical differentiation 
we computed |uo), which is solution of the eigenvalue problem for k = 0. 

We also note that the autocatalyst exhibits a slowly decaying tail whose length 
increases as k approaches zero. The length of the tail can be estimated from the equa- 
tion, X2 = k'^X2 — c(l — /w)~^X2) "which is obtained from the equation governing the 
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dynamics of the autocatalyst (second member of by setting the reaction terms 
to zero. When k is small the solution exponentially decays as exp 

(-(1 -/i)A;Vc) 

and the length of the tail is of order Lta ~ c[(l — 
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FIGURES 

FIG. 1. Panels (a) and (b) are space-time plots of the minima of the autocatalyst (B) 
profiles for the diffusion coefficient ratios 5 = 5 and 6 = 8, respectively. The isoconcentra- 
tion level for these plots is l3r = 0.5. The ordinates range from to L, where the system 
size is L = 2048. In the abscissas the initial time is t = and the final time is t = 100 000 
for (a) and t = 200 000 for (b). 

FIG. 2. The figure shows C{t,y) for 5 = 5 for several values of t indicated on the 
plot, as well as C{y) and the theoretical value of C{y) = Ty{L — y)/T>L for a diffusive in- 
terface where nonlinear terms in the Kardar-Parisi-Zhang equation can be neglected. The 
parameter T) = T /T), determined from a fit to the C{y) data, has the value T) « 0.2. The 
parameter F may be estimated independently from the linear growth of the A: = Fourier 
mode. Its value was found to be F ~ 0.018. Fronts corresponding to the isoconcentration 
lines [3r = 0.5 were used to construct these correlation functions. The system length was 
L = 256 and C{t,y) was determined from an averages over 25 realizations of the front 
evolution. In the computation of C{y), we used to = 25 000 and to + T = 50 000. 

FIG. 3. Power spectra for 5 = 5 and 5 = 8. The isoconcentration lines for (3r = 0.5 
were used to construct these power spectra. The system length was L = 2048 and E{k) 
was determined from a single run of the front evolution. For 5 = 5 we used to = 25 000 
and to + T = 100 000, while for 5 = 8 we used to = 50 000 and to + T = 200 000. 

FIG. 4. Space-time plot of the minima of the profiles for (a) the Kuramoto-Sivashinsky 



equation and (b) the amplitude equation (19) 



FIG. 5. (a) The solid line is a schematic representation of the dispersion relation for 
5 > 5c- The dotted line shows the hypothetical form of the dispersion relation with two 
energy sources, (b) Proposed rate of energy mixing g(k) that produces biscale patterns. 
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FIG. 6. (a) Plot of the potential in the WKB approximation, (b) Plots of the 
components of the left eigenvector for the case of equal diffusion coefficients. 

FIG. 7. Dispersion relation of the cubic autocatalysis model for the case n = 0.45 
(or 6 = 2.636). The symbols '+' denote the dispersion relation constructed from the 
numerically obtained growth rate. It is compared with the first order approximation 
uj{k) 0.142A;^ for small values of k. 
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